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HIGHER-ORDER FINITE ELEMENT METHODS FOR ELLIPTIC 
PROBLEMS WITH INTERFACES 

JOHNNY GUZMANS MANUEL SANCHEZ-URIBE^ AND MARCUS SARKIS ^ 


Abstract. We present higher-order piecewise continuous finite element methods for solving a 
class of interface problems in two dimensions. The method is based on correction terms added 
to the right-hand side in the standard variational formulation of the problem. We prove optimal 
error estimates of the methods on general quasi-uniform and shape regular meshes in maximum 
norms. In addition, we apply the method to a Stokes interface problem, adding correction terms 
for the velocity and the pressure, obtaining optimal convergence results. 

Keywords-. Interface problems, finite elements, pointwise estimates. 


1. Introduction 


In this paper we continue the work started in |24j and consider higher-order piecewise continu¬ 
ous finite element approximations to the following interface problem; Let C be a polygonal 
domain with an immersed smooth, closed interface T such that U = U U n"*" and T encloses 
Consider the problem 


(1.1a) 

-Au = f 

in 

(1.1b) 

u = 0 

on dit, 

(1.1c) 

[u] = 0 

on T, 

(l.ld) 

[DnU] = 13 

on T, 


where the jumps across the interface F are defined as 

[u] = u'^ — u~ , [Dnu] = Dn-u~ -|- Dn-eu^ = Vu“ • n~ -|- Vu'^ ■ . 

Here we denote hy = u\q± and is the unit outward pointing normal to 

Numerically, the problem is to find an approximate solution on meshes not aligned with the 
interface, that is, we allow the interface to cut elements. In this context, the finite difference 
methods by Peskin [301 [31] (i.e. immersed boundary method) and by LeVeque and Li [27j (i.e. 
immersed interface method) are the most renowned. Both methods were developed for more 
involved problems and for lower-order finite differences techniques. The aim of this paper is to 
develop higher-order methods based on finite element methods and to establish a priori pointwise 
error estimates. We consider the Poisson interface problem 0 and the Stokes interface problem 


(1.2) which will be fundamental towards developing very accurate methods with a rigorous finite 


element analysis for more involved problems. 

Naturally, finite element versions of the methods above have appeared; see for example 
[I8lia[ll|2llini[8lllll[l2l|9lll3l[2ll|23l[22l|2l!. In our recent work |24j we derived a piecewise 
linear finite element method for the above problem and proved it is second-order accurate. The 
attractive feature of the method in |24j is that only the right-hand side needs to be modified, 
which is one of the advantages also of the immersed boundary method (see m,m) and immersed 
interface method (see [26], |3|, [20) 1 for the above problem. Moreover, the correction term added 
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in 


is only based on the edges intersecting the interface, this is due to the fact that test 
functions, linear polynomials, are harmonic. However, it seems that the approach used in [23] 
cannot be generalized to higher order approximations. Therefore, in this paper the approach we 
take is based on correction functions that live on the entire triangle instead of just correcting 
terms on the edges. 

Motivated by fluid applications, we naturally seek methods for higher-order finite element 
spaces, and in this context edge based modifications are not enough. Guided by our recent work 
in the piecewise linear case and also by other papers; see for example [20l|T8], it appears to us 
that the construction and addition to the right-hand side of a correction function is the key to 
achieve a higher-order method. This also appears in the finite difference context, for example 
a fourth order method was developed by Marques, Nave and Rosales [29] using a correction 
function approach. Very recently this idea was materialized by Adjerid, Ben-Romdhane and 
Lin [T]. They developed higher-order methods for problems involving discontinuous coefficients 
(which is a more general problem). However, they use strongly the assumption that the interface 
is a straight line. The key is to use both, the jump condition and the PDE, to find higher-order 
jump conditions. Inspired by their results, we define correction functions for any polynomial of 
degree k for curved interfaces. 

The contribution of this paper is in the direction of [I], we develop a higher-order piecewise 
continuous finite element method for problem (1.1) with curved interfaces. Specifically, in Sec¬ 


tion]^ we develop notation and propose a finite element method, for each polynomial of degree 
k, introducing a correction function and adding it to the variational formulation, only mod¬ 
ifying the right-hand side of the equation. To do this, we construct this correction function 
incorporating the jump conditions of the exact solution on the interface. To the best of our 
knowledge, this is the first family of numerical methods (with any order approximation) were 
one can prove optimal accuracy for the above problem without modifying the stiffness matrix. 
Besides the novel method, an important contribution of this paper is the techniques used to 
give a-priori error estimation analysis for the Poisson and the Stokes interface problems. In 
particular, an interpolation estimate using the correction function is given(see Lemma]^, which 
is crucial for the full analysis of the proposed method. In Section]^ we prove that our method is 
k + 1 (mod a logarithmic factor) order accurate in the maximum norm, if piecewise polynomials 
of degree k are used. 

As mentioned before, in this paper we also consider a finite element approximation to a Stokes 


interface problem, i.e., under the same geometry assumptions for problem (1.1), we seek for a 
velocity vector u and pressure p satisfying 


(1.2a) 

—Au + Vp = / 

in Q, 

(1.2b) 

V ■ u = 0 

in n. 

(1.2c) 

u = 0 

on d^l 

(1.2d) 

[BnU -pn]= (3 

on P, 


As is well known, time dependent versions of problem (1.2) have many applications in biology. 


see for example M- In fact, the immersed boundary method is used primarily for solving time 
dependent versions of (1.2) with possibly nonlinear terms; see for example Peskin and Tu 


Subsequently, LeVeque and Li m consider this problem applying immersed interface method 
techniques. Some extension and finite element versions of these approaches can be found in 


|5l|25l|28|. We consider the numerical method in this paper for problem ( |1.2[ ) as an important 
step toward defining higher-order methods for time-dependent problems with moving interfaces. 
We will show in Section that the same methodology used to achieve higher-order methods for 
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Poisson problem (1.1) can be applied to Stokes problem ( |1.2[ ), achieving optimal convergence 
results. 

In Section we test the methods introduced in Section and with some numerical examples 
that illustrate the properties proven in previous sections. Furthermore, we provide in the Appen¬ 
dix quadrature formulas for the integration over curved regions crucial to achieve higher-order 
results, although our analysis considers exact integration. 


2. Finite element methods 


In this section we present a finite element method for problem (1.1) using continuous piecewise 
polynomials of degree k. We assume that the data jS is smooth. Furthermore, we assume that 


u 


± G and /b± = G 


2.1. Notation. Let Th, 0 < /i < 1 be a sequence of triangulations of II, 12 = UreT^P, with 
the elements T mutually disjoint. We assume the mesh is shape regular, see [7]. We adopt the 
convention that edges, elements, regions are open sets, and we use the overline symbol to refer 
to their closure. Let denote the diameter of the element T and h = maxT’/i'p. Let Vh be the 
space of continuous, piecewise polynomials of degree k, i.e., 

Vh = {v^ C(Q) n B^(Q) : vjr G F^(T) VT G Th}, 

where P^(T) is the space of polynomial of degree less than or equal to k on T and IIq(Q) the 
space of functions in vanishing at dfl. 

Next, we define an interpolant onto Vh- 

Definition 2.1. Given G we define locally Tv G Vh sueh that 

GD = {lS: 

for all 6 gT, the degree k Lagrange points ofT. 


Note that if v is continuous Tv is simply the Lagrange interpolant of v. However, if v is 
discontinuous then Tv interpolates values of v on Lagrange points not intersecting F and for 
Lagrange points lying on F it takes the values of v coming from Ll~ (this is without loss of 
generality). The following proposition states the stability result of the interpolant T- 

Proposition 1. Let G C'(II^) and Ih defined above, then we have 
(2.2) < C\\v\\Loo(^rp'j VT G Th- 

To prove Lfi based error estimates shape regularity of the meshes suffices. However, as is well 
known, some form of quasi-uniformity of the mesh is needed to prove max-norm estimates, even 
if there is no interface [?]. Since we will prove max-norm estimates of our method, and in order 
to avoid unnecessary details, we will assume that the meshes are quasi-uniform. 

Also, for the sake of making the presentation simpler to the reader, we make the following 
assumption. 


Hypothesis 1. We assume here that the interface F intersects the boundary of each triangle 
TGTh at most at two points- If F intersects the boundary of a triangle T in exactly two points, 
then these two points must be on different edges e ofT. 

Next, let 7)f denote the set of triangles T G Th such that T intersects F, that is, T n F / 0. 
For each T G T^ let yT and zt be the two endpoints of T n F, see figure Let Lt denote the 
line segment connecting yT and zt- Let ij be the unit vector perpendicular to Lt and pointing 
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outward T . Let also r be the unit vector parallel to the line Lt such that r is the rotation of 
rj by ninety degrees counterclockwise. 

For each i = 0,1,... ,k, let denote the Gauss points of the segment Lt- For each 

X-’ , let Q-’ be the line perpendicular to the line segment Lt that passes through the points 
X-’ . We then define x J n F, for i = 0, ..., Note here that the choice of Gauss points 

is a preference of the authors, related to the quadrature rules, but not essential in the proofs. 
We could also use, for instance, equally spaced points. 



Figure 1. Illustration of our notation, for T G 7^^. 

Letting := T n we define the following space for T G 
(2.3) S^{T) = ju; G L^{T) : w\t± G P^(r^)} . 

2.2. The proposed finite element method. We now present our finite element method for 


problem (1.1). Find u/j G 14, such that 
(2.4) 


Vu/i • Vu dx = 


[ fvdx+ / jdvds — 

n Jr 


E 

rerp' 


S/w^ ■ Vvdx 


for all u G 14 , where is a correction function to be defined in Section 2.3 


2.3. The correction function. We now show how to construct a piecewise polynomial function 
G S^{T) that will help to correct the right-hand side of the natural finite element method 
((2.4) without the correction term) to render it higher-order. We note that the functions w^, 
for T G , are discontinuous across elements and satisfy jump conditions at Gauss points on 
F n T. Suppose that we give you a function u, then we let be the unique function in S^{T) 
(see Lemma that satisfies 


d!^ ^Wt{x\'^) = ) for 0 < i < f and 0 < I <k, 

Wt{6) =0 for all degree k Lagrange points 6 of T. 




(2.5) 

( 2 . 6 ) 
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We would like to stress that if u is the solution to (1.1) then we know all the jumps of u in 


terms of the data j3 and / (see Section 2.4) so we can construct a priori. It is important to 
notice that we impose directional derivatives jump conditions in the 77 -direction rather than in 
the n-direction. A reason for this choice is because we can show unisolvence, when F is curved, 
for any polynomial of degree k. Additionally, the construction of the correction functions can be 
made explicitly, that is, it does not require solving any local linear system. Such construction 
permits us to develop not only an elegant analysis for any polynomial degree k but also an 
analysis without imposing strong condition on how F intersects an element T in order to control 
the ill conditioning of this matrix. We next show how to transform jump data described in 
n-direction to jump data in the 77 -direction. 


2.4. Data of jumps. In this section we show that we can obtain the jumps of higher derivatives 
of u across F from data, (3 and /. Similar ideas were used in [2 ^ 1 ^ 1^. 

Fix a point x G F. Let n be the normal vector to F at x, and t the tangent vector to F at x. 
We clearly have [Dnu{x)] = /3 and [Dtu{x)\ = 0. In fact, we also have [Zl|u(x)] = 0 for any 
Note that from the Poisson’s equation (1.1a), we have 


-D^u- Dju = f inO=^, 

and so [D^u(x)] = [/(x)]. Moreover, we note that [DtDnu{x)] = Dt [Dnu{x)] = /3'(x) on F. 

Now we proceed by induction. Suppose we have all jumps of the derivatives of order I — \ 
in terms of f3 and /. Then, we will show how to get the jumps of derivatives of order L Let 
i + j = i- If 7 > 1 then 

[DlDiuix)] = Dt [Dl^Diuix)] 

Otherwise, using Laplace’s equation we have 


D^u{x) = f{x) 


-Dt 


DtD^-^u{x) 


Suppose that we would like the jump of tt in a different direction, say 77 . Then, we write 
77 = an + bt obtaining 


d1^u{x) 


1=0 


D^dI ^n(x) 


3. Error Analysis 

The objective of this section is to prove rigorous pointwise error estimates for the above 
method, which we achieve in Theorems andBefore doing this we need some technical results 
associated to the subspace S^{T) and approximation properties of the correction function wif,. 

3.1. Properties of S^{T). We now introduce some crucial lemmas related to the space S^{T). 
Lemma 1. Let T G and define tt = \Lt\- For any v G P^(T), we have 

k 

max |D^"^u(xf'^)|, j = 0,1, 

e=o 

where the constant C depends only on the shape regularity of T, the polynomial degree k and the 
regularity o/F. 
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Proof. Let be the segment with length \L^\ = 2hT, centered at {i/t + zt)/‘2 and aligned with 
Lt- This guaranties that any point x G T can be projected orthogonally on L^. Then, using 
Taylor’s expansion on T from L^, we can easily show that 

k 

h^r\\D^vh^(^T) < CY,4~'\\D’:,-^v\\L^iL!fy 
e=o 

Using Taylor’s theorem on from a point of Lt we get 


j = 0,1. 

The inverse inequality gives 

\\D%D'f^-^v\ 

We therefore have 


s=0 


(j 


4II 






1=0 s=0 

. 1 1 


k—l II 

’r) ^IIl°°(Lt) 


c'4EEi?LLiiob'>’iii»(ix) 

e=o s=o '^T ’’t 


k I 


s c4EE;i^ll^L'”IU"(iT) 

i=0 s=0 T 

k 


c^Eyiia 


e=o 


k—i II 

T) ^IIl”(Lt)’ 


where we used that tt < hT- 
To bo 
we have 


To bound the right-hand side above, we use induction on i. First, using that D^v is a constant 


|Zl>||i^(i^) = |Il4(xE)|. 


Assume that we have proved 

k ^ 

( 3 - 1 ) E 3 


1 


- II T^h—f II 

^ rL ^ rf o<i<i' 

t=m+l ^ (.=m+l ^ 


< C max \D^ ^v(x 




then we want to prove that 

k k 

(3.2) E ^WD’^r^vh^iLr) < C ^ ^ m^^\D>f^v{xf)\. 

i=m ^ £=m ^ 

Since D^~^v is a polynomial of degree m, we have that 

l|Ob”«lll»(lrt < C 

' 0 <^<m ' 

Using Taylor’s theorem we have 

k 

\L>r, v{x- }\<\v^ u(x. j|+ 2^ « ||L»^ X||Loo(ij,), 

i=m+l 



HIGHER-ORDER FEM FOR ELLIPTIC PROBLEMS WITH INTERFACES 


7 


where d = max It is clear that d < Cr^, since we have assumed that F is 

o<i<e<k ' * « ' 

smooth, that is the radius of curvature is 0 ( 1 ), however we note that we will only use that 
d < CrT- Hence, we have 

k 


^ max ^ max 

r™ 0<i<m ' * r™ 0<i<m ' 


' rjn 

However 


■o(if’^)|+C ^ 


1 




C=m+1 T 


1 „2m 
T 


^ for £ > m + 1 and so using (3.1) we arrive at (3.2). 


□ 


The following is a fundamental lemma for the construction of and the estimation of 
hip\\D ^when we choose 




, see equation (2.5). 


Lemma 2. Given data {ci^i} for 0 < i < £ and 0 < £ < k. There exists a unique funetion 
w G S^{T), such that 


(3.3) 

(3.4) 


*^w{xf'^) = Ci^i for 0 < i < £ and 0 < £ < k, 

w{6) = 0 for all the degree k Lagrange points 0 ofT, 


with the following bound 

(3.5) < 


^ ^ 1 

Y] -f lo,£|, for j = 0 , 1 , 
rL o<i<£ 

1=0 T — 


where C depends only on the shape regularity ofTh, the polynomial degree k and the regularity 
ofV. 

Proof. We will construct w of the form w = z — IhZ, where z G S^{T). Notice that, by definition 
of Ih, z — IhZ vanishes on the {k + 1)(A: + 2)/2 Lagrange points of T, satisfying (3.4). Moreover, 
since IhZ is smooth on T 


k-l„ 


u;(xf^) = 

1 

Q 

1_1 

Jo 

in r+ 

z = < 



in T~ 


i / ei\ 
z{xf ) 


D: 


The function z will be given by 


where v is the unique polynomial on P^, such that 

ioi £) <i < £ and £) < £ <k. 

The existence and uniqueness of v follow from representing v = v{s,r) as a polynomial of 
degree k in r (r-direction) and s (ry-direction), where s = 0 represents the straight line passing 
through yx and zt, then by decomposing 

v{s,r) =Pk{r) + spk-i{r) + s^pk- 2 {r) H-h s^o, 

where for £) < £ < k, is a, polynomial of degree £ in r. It is easy to see, by using interpolation 
—0 T 

at the Gauss point Xg’ , that po exists and is unique, then by using interpolation at the Gauss 
points x^*’^ and x\^ that pi(r) exists and is unique, and so on. 

According to Lemma we have the following bound 

^ ^ 1 

(3.6) hk.\\D^ v\\x°°ix) < C —f max \ciA. 

' r^ o<i<l ’ 
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Hence, w satisfies (3.3). Moreover, 


“ hz)\\L°°{T±) < ^Tll^^^llL“(r-) + ^tII^^42;||l°=(T)- 

Using an inverse estimate and stability of the interpolant, we have 

^t\\D^I hz\\L’^(T) < C\\Ihz\\L<x(^j^^ < C ||2;||ioo(7’) < C ||r;4oo('r-). 

Hence, 

U IIIIl^(T“)• 

We get (3.5) once we apply (3.6). 


□ 


3.2. Approximation properties of w^. Since we are assuming G and T is 

smooth, there exist extensions G (see Lemma 6.27 [IS])> such that the following 

holds 


U 77 - U 


± 


onH=^, ||u|||cJ=+i(n) < C'||u=^||cfc+i(n±) 


±1 


Let i? 2 rj, C r® be a ball of radius 2rT that encloses . Here is the smallest patch of 
triangles of the mesh Th on the neighborhood of T G . Let Jt be the projection onto 

polynomials of degree k in H 2 rj, and consider its natural extension to all of . Then, we can 
prove the following lemma. 

Lemma 3. Let w G then we have 

(3.7) 4||D4Jt(iu) - 'w)IU-(B2,^) < C r^'^^\\w\\ck+^B2rj,) 0 < j < k, 

and 

(3.8) hip\\D^ {Jt{w) — r(;)||£oo(2’S) < C /i^^^||iu||c*+i(T^) 0 < j < k. 


Proof. The inequality (3.7) is a standard approximation of the projection. To prove (3.8) we 
apply Taylor’s theorem to get 


k-j 


h^j.\\D^{J t{w) - w)\\l<^(^te) < h^T'^h^B\\D^^\jTiw) - w)\\L^^B2rj.) + 


i=i 


The result follows after applying (3.7) and using that vt < hr- 


□ 


The following lemma establishes the approximation result for the correction function defined 
in Section! 


Lemma 4. Suppose the solution u to problem (1.1) satisfies G C'^'''^(H^), and wf is the 
correction function defined by (3.3) and (3.4). Then, we have 


^tII-^^^ “ — 'u^T)llL°°(r±) ^ 4 hjf^ (^II^Bllc''+i('r®) ll'^^ullc'^+Tr®)j > i — Oj 1 

and C depends only on the shape regularity of T, the polynomial degree k and the regularity of 

T. 

Proof. We will define v G S'^(r) (see ( |2.3[ )) as follows 

f Jt{ub) on T+, 

|Jr(tt^) onT”. 


V = 
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Clearly v — I^v = by Lemma Hence, we have 
h{,\\D^{u - Ihu - wt)IIl“(t±) = h{.\\D^ ((n - u) + (t> - hv) + Ih{v - u) - w^) \\l^(t±) 

= h{.\\D^ {{u -v) - Ih{u -v) - ||z,°°(r±) 

< C {\\D^{u - i;)hoc(r±) + \\D^wlf^hoo^r±)) 

+ C \\u — w||Ltx>(2’), 

where we used — w^, an inverse estimate and the stability of Ih in the max-norm. 

Using (3.8) we get 

h^T\\D^{u — -|- IIu — u||x,oo(7’±) < Ch’p~^{\\u^\\Qk+i(^rpE'^ + ||w^||c'=+i('r®))' 

Estimate ( |3.5| ) implies 

k 


h^\\wlf^\\L^^T) <Chl^Y, 


1 


e=o 


k-e I 


D^iv - u) 


lL°°(rnr)' 


Applying (3.7) we obtain 

k 


E 

e=o 


1 


,k-e I 


dI^{v — u) ||i,<=o(rnr) < CrT{\\u'^\\(jk+i(^B^ \ + \\u^\\ck+T-(B2rj,))’ 


which completes the proof. 


□ 


3.3. Error estimates. The next lemma will show that the correction term in the finite 
element method (2.4) will allow us to compare I^u — u^- 


Lemma 5. Let G be the solution of (1.1) and he the solution of (2.4). Then, 

it holds 

^ V(4t6 - Uh) • Vudx < Ch'"\\Vv\\Li(n) (lk'^llc'=+i(f^+) + ll^~llc'=+i(f^-)) ^ 

where C depends only on the shape regularity of {Th}h>o, the polynomial degree k and the 
regularity ofV. 


Proof. 


V{IhU — Uh) ■'^v dx = / V{IhU — u) ■ Vv dx + / Vu-Vvdx— / Vuh-Vvdx 
! Jn Jn Jn 

= / V{IhU — u) ■ Vv dx + / Vwlf-Vvdx 


T&Tf ■ 


= / V{IhU — u)-Vv dx + / V{u — IhU — wlf) ■ Vv dx 

Ten\Tf Terr 

The result now easily follows from Lemmaand the fact that u is smooth on T G Th\T^■ □ 

From the above lemma we can easily prove an optimal estimate in the semi-norm: 

||V(4m — Uh)\\L2(n) <Ch^ (ll^'''llc''+i(f2+) 11'“ llc''+i(f^“)) • 

However, are goal is to prove estimates in the maximum-norm as our next result states. A 
slightly sub-optimal (off by a log factor) can be proved if we use the above lemma directly. In 
order to prove the optimal estimate, we will give a more involved argument. 
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Theorem 1. Suppose that O is convex. Let G be the solution of 

the solution of ( |2.4| ), then 

\\V{IhU - u/i)||loo(q) < C /i^(||ii^||c"=+i(n+) + ll'“~llc'=+i(f^-))’ 

where C depends only on the shape regularity of {Th}h>o> the polynomial degree k and the 
regularity ofT. 

Proof. Let eh = IhU — Uh and suppose that the maximum of IdxiChl occurs at z G (for some 
fixed 1 < i < 2). Suppose z £ T, for some T G Th- Consider now the regularized Dirac delta 
function = 5h £ Cq{Tz) (see m), which satisfies 

(3.9) r{z) = ir,6h)T., W £ P'^{T,), 
and has the following property 

(3.10) \\hh\\w-:i(T,) < 1 < g < oo, r = 0,1. 

For each i = 1,2, define the approximate Green’s function g £ which solves the 

following equation: 


(1.1) and Uh be 


(3.11a) -Ag = d^iSh in D, 

(3.11b) g = 0 on dLl. 

We also consider its finite element approximation gh £ Vh that satisfies 


/ ^9h ■ Vu dx = V ^x^6h dx for all v £ Vh- 

In Jn 


(3.12) 

From the work of Scott and Rannacher |32j we have 

(3.13) \N{g-9h)\\LHn)<C. 

Moreover, using a dyadic decomposition one can show 

(3.14) \\Vgh^n)<Clogil/h). 

A log free estimate holds if we consider a smaller domain, i.e. 

(3.15) l|V5ki(sr) < C, 

where S''" = {x G D ; dist{x,r) < uh} for some fixed constant k; see for instance 
combining (3.13) and ( 3.15| ) we have 

(3.16) l|Vfl'/i||£,i(5r) < C. 


We start by using the dehnition of 6h and problem (3.12) 


II||_L°°(o) — l^xi(. 2 )I — I / 6hdxi^hdx\ — I / Ox^dhehdx\. 

Jn J n 


Then, we see that 


\\dxieh\\L°°(n) = \ / V9-Vehdx\ = \ / Vgh -Vehdx]. 


. Hence, 


If we follow the proof of Lemma we see that 
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where 


and 


Ji = ^ / V{IhU- u) -Vghdx, 


TeTh\T^ 


J2 = ^ y V{IhU +- u)-Vghdx. 


Terr 


Applying Cauchy-Schwarz inequality to J 2 we get 

J2 < C\\Vgh\\Li{s^) Toxax\\V{IhU + WT-u)\\L^(^T)- 




Moreover, using (3.16) and Lemma we have 

J2 < C h^{\\u^\\(jk+Hyi^+-^ + ||tt 

To give an estimate for Ji we define Rh = U^^gT-^^T-rT. Now, adding and subtracting V 5 , we 
obtain 

Ji = ^ V{IhU-u)-Vghdx= {V{lhu-u)-Vigh - g)+V{IhU-u) ■Vg)dx. 

Jt Jr. 


reTArr 


Using (3.13), we have 

/ V{IhU-u)-Vigh-g)dx < Ch^{\\u-^\\ck+i^Q+) + \\u-\\ck+i^Q-)). 

J Rh 

For the remaining term we integrate by parts to get 

/ V{IhU-u)-Vgdx = / {IhU- u)dxi5hdx + / {hu - u)Dng ds. 

JRh JRh JdRh\an 

Here we used that, since H is convex, V 5 is continuous and so integration by parts makes 
sense. 

Clearly we have 


' Rh 


{IhU — u)dxiShdx < Ch^ ^Il^^llc'=+i(r2+) + 11^ Ilc'=+i(r2-)) • 

Finally, we have 

{IhU - u)Dngds <C\\IhU - u\\L^(^ii^)\\Dng\\L^gji^\QQ) 

<Ch^~^^\\Dng\\Lil^gii^\gQ) ||c*+i(n+) + \\R llc^+qn-)) 


ldRh\dn 


In Appendix we prove the bound 


C 


\Dng\\mdRh\an) < 


(3.17) 

which will then show that 

Jl < Ch’^{\\u'^\\(gk+l^Q+-j + ||tt ||(^fe+l(f^-)), 

and will complete the proof. 


□ 


Next, we will prove an estimate for the error in the maximum norm. In the case there is no 
interface a logarithmic factor is not present for k > 2 (see 0 ), however, we do not see how to 
remove this factor in our setting. 
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Theorem 2. Suppose that O is convex. Let G be the solution of 

the solution of ( |2.4| ), then 

(3.18) \\IhU - Uh\\L^(n) < Ch^^^ log{l/h) (^||n+||c.;=+i(j^+) + ||m" ||cj=+i(q-)) • 

where C depends only on the shape regularity of {Th}h>o, the polynomial degree k and the 
regularity ofT. 

Proof. We follow the proof of Theorem 5 in [2^ . Let z ^ LI he arbitrary and let 5^ = defined 
in proof of Theorem Let g satisfy 


(1.1) and Uh be 


(3.19) —= 5h in LI, 

(3.20) g = h on dLl, 

and consider its continuous piecewise linear finite element approximation g^. Then 
{IhU-Uh){z)= / {IhU- Uh)Shdx = / V{IhU-Uh) -Vghdx 


T&n\rH 

= : Ji + J 2 

We first give a bound for J 2 


V] / V{IhU-u)-Vgh)dx+ / V{IhU + wf-u)-Vgh)dx 

Jt JT 


Terr ■ 


J 2 < C m.a:!c\\IhU + wf-u\\Loorrp)\\Vgh\\LHT) 

T&Tf 

< Ch^ (ll^‘’''llc'=+i(!^+) + 11“ Ilc'=+ 1 (!^")) IY5/i|lLi(sr) 

< Ch’^ (^ll^^’'’llc'=+i(!^+) + 11“ Ilc'=+i(f1")) + ||V5||j;^i(5r)^ . 

In [23] we proved < Chlog{l/h), therefore 


•h < log(l//i) (^||u^||c'=+i(r2+) + 11^ Ilc'=+i(r2-)) • 

Now for Ji, we will use the Raviart-Thomas projection (see |33]) LI : H^{Ll) —)> defined 

locally for any T G 7/t, nlT : H^{T) —)• RTq{T), where 

RTq{T) = [P°(r)]2 © xP°(r), 4 .^ = {(^ G [L^{Ll)f : ct)\T G RTo{T) VT G Th] . 

Then, we observe 

Ji = Ji{Vgh) = JiiVgh - LlVg) + Ji{UVg). 

In the Appendix of [21] we prove the estimate \\Vgh — nV 5 ||j;,i(Q) < hlog 1/h, then we clearly 
we have 

>^i(V5h —nv^) < Ch’^ (^\\u'^\\Qk+i(^Q+'j + \\u Ilcfc+Lrj-)) “ nVg||j;^i(Q) 

< log(l//i) ^||u’''||( 3 .fc+i(Q+) + ||tt . 

Using that n(V 5 ) is piecewise constant and has continuous normal components across edges, 
we have after integration by parts 

Ji(nV5) = ^ l'{IhU-u)UVg-n, 

e 
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where are set of edges that are both an edge of a triangle in Th\T^ and a triangle in ThXT^ 
Therefore, we see that 


L 


Ji{YiVg) = ^ hlhU-u)IWg-n 


eGS, 


r,d ' 


< Ch^ ||nV5||j;^i(5r) ^||w’'"||c'=+i(D+) + 11“ Ilc'=+1(!^")) 

<Clog{l/h)h^^^ ^||tt"''||(^fe+i(f^+) + ||tt , 


where again we used ||nV 5 ||ii( 5 r) < Chlog{l/h), which follows from results in [2l] . 

□ 

Remark 1. Since u — {uh + = {IhU — Uh) + (u — — Ihu), by using triangle inequality and 

Lemma^ Uh + wl^ approximates u on T by the same estimates given in Theorems^ and[^ 


4. Stokes Interface Problem 

In this section we consider the Stokes interface problem in two dimensions introduced in 
equation (1.2). Equivalently, we can incorporate the jump condition, equation (1.2d), as follows 


(4.1a) 

(4.1b) 

(4.1c) 


—Alt + Vp = f + B in R 
V • It = 0 in II 
It = 0 on dQ 


where 


B{x) = f P{s)6{x — X{s))ds, 

Jo 


and Ai(s) with 0 < s < ^ is the arc-length parametrization of the interface P. 

Following the ideas of LeVeque and Li in 122!, we can easily write individual jump conditions 
for the velocity and the pressure in terms of the tangential and normal component of the data 
(3. Let 9 be the angle between the x-direction (x-axis) and n-direction pointing outward the 
interface P at a point X{s). Then, we write the normal and tangential components 


Pis) = 


Pi 


cos(0) sin(0) 


Pis). 


— sin(0) cos(0) 

The jumps conditions for the velocity and the pressure are given by 


(4.2) 


b] = Pu 

[DnP] = 


[li] = 0, 

[i:»„,ii] = P + Pin. 

For the sake of completeness, we present the derivation of this jumps in Appendix [B) 
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4.1. The finite element method. 

Stokes interface problem (1.2). Find 


We first present the standard variational formulation of 
'u,p) G [Hq{Q)]‘^ X Lq{Q,), such that 


(4.3) / Vu-.Vvdx— / pV ■ vdx = / f ■ vdx + / jS ■ vds Vu G [Ffg(n)]^, 

Jo. Jn Jn Jr 

/ qV ■ udx = 0 V(? G Lg(n), 


where = {q G L^(n) : f^q = 0}. 

As before, we consider a sequence of triangulations of Cl, Th with 0 < /i < 1, and 11 = 
with the elements T mutually disjoint. Let /it denote the diameter of the element T and 
h = max'T Ht- We assume that the mesh is quasi-uniform and shape regular. 

We consider a class of finite element subspaces Vh C [iLQ(ll)]^ and Mh C Lq{Q,) satisfying 
the following assumptions: 

A1 Vh and Mh are a pair of inf-sup stables subspaces, with Vh C [F1q(11)]^. 

A2 We let /c > 1 as the maximum integer such that 


Vi := 


:= [v G C(ll) n : ^It G VT G %} C Vh, 

and, if Mh contains the discontinuous pressure space of degree A: — 1 we let 


M; 


otherwise 


:= [q G Llin) : q^ G VT G %} C Mh, 


:= {q G C{n) n Llin) : q\T G VT G %} C Mh- 


M, 


For instance. A: = 1 for the pair ~ reduced ~ and mini element, while k = 2 for 
Taylor-Hood P| — Pi. 

We next define the interpolant onto these spaces. We let Ih be the interpolant defined corn- 


space M^ in the case of continuous pressure finite element spaces. Otherwise, if Mh contains 
discontinuous finite element pressure spaces we define Jh to be the projection onto M^~^. 


ponentwise in (2.1) onto the space V^, and Jh be the interpolant defined in (2.1) onto the 


Find {uh,Ph) G Vh X Mh, such that 


/ Vuh'.Vvdx— / phV ■ vdx 
/ n Jn 


/ / • vdx + /3 ■ vds 

In Jr 


(4.4) 


T&TJ 


: Vu dx+ • vdxj , 


[ qV-Uhdx =- V / qV-w^dx, 
Jgi JrjJr Jt 


rerp ■ 


for all {v,q) G Vh X Mh- 

Using the last two equations of (4.2), the correction function is defined componentwise as 
in Section 2.3 The same is true for when Jh is the Lagrange interpolant by using the first 
two equations of (2.3). In the case Jh is the projection we replace equation (2.6) with the 
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condition Jh{w^) = 0, 


i.e. 


(4.5) 

(4.6) 






Pi^i 


ior 0 < i < £ and 0 < i < k, 


Jh{w^) =0. 


Note that each component of and are obtained independently to each other, therefore, 
the Lemma can be applied to each one separately. 

Similarly to Lemma we have the following estimate for the interpolant Jh and the correction 
function 

\\p — JhP — < C ^l|p^llc'=+i(r) + 11^* llc''+i('r)) ^ '^h ) 

where C is a constant depending on the shape regularity of T, the constant k and the regularity 
of r. Then, the following result holds and can be proved similar to Lemma 

Lemma 6. Let {u,p) be solution of (1.2) and assume that G [C^'’“^(n^)]^ andp^ G C^{Ll^). 
Let Vh and Mh be the finite element spaces satisfying assumptions A1-A2 and consider the 
definitions above for k, Ih and Jh- Let {uh,Ph) G Vh x Mh be solution of ( |4.4[ ). Then, it holds 

/ V{IhU-Uh) :Vvdx - / {JhP - Ph)V ■ vdx 

J V. Jo. 

< C/i^||Vr^||2,i(o) (^||u’''||c'fe+i(D+) + 11'“ Ilc'“+i(t2") lll^^llc''=(o+) + Up IIc''(d-) 


In 


qV ■ {uh - Ihu)dx < Ch^\\q\\Li(n)(\\u-^\\ck+i(n+) + \\u llcfe+Tf^-)) ’ 


for all {v, q) G x Mh- 


Proof. For the first equation we use (4.3) and (4.4) to obtain 


in 


V{IhU-Uh) :Vv dx - / {JhP - Ph)V ■ vdx 

Jn 

= / V{IhU — u) : Vv dx — / {Jhp — p)V-vdx 
Jn Jn 

+ ^ ■ vdx + / Vwlf:Vvdx] 

Terf ^ 

V{IhU — u) : Vv dx — J {JhP — p)V ■ vdx^ 


T&Th\Tf 


Terf 


+ ( / V {IhU — u + wlf) : Vv dx+ {JhP — p + w^)V ■ vdx j . 


Similarly for the second equation we have 
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Then the results follow from the properties of the correction functions and Lemma 0 


and the fact that (u, p) is smooth on T G Th\Tl 




□ 


Analogously to the proof of Theorem we can prove the following result using approximate 
Green’s function estimates for the Stokes problem. We give a sketch of the proof. 




G 


Theorem 3. Let {u,p) be solution of (0 and assume that G and p 

C^(n^). Let Vh and M/j be the finite element spaces satisfying assumptions A1-A2 and consider 


the definitions above for k, Ih and Jh- Let {uh,ph) G x be solution of (4.4). Then, there 
exists a constant C > 0, such that 




(4.7) ||V(//iW — + II J/i(p) — < Ch^ (^\\u'^\\(jk+i(^Q+'f + \\u 

+ l|P^llc'=(o+) + \\P \\c>={n-) 

Proof. Let = IhU — Uh and suppose that the maximum of dxj{eh)i {i,j = 1,2) occurs at 
z G n. Suppose z £ T £ Th- Considering the definition and properties of the regularized 


Dirac delta function 5h = df introduced in (3.9) and (3.10), we define the approximate Green’s 
function {g, A) G [L7q(D)]^ x L‘^{LI) such that 


—Afif + VA = {dxj6h)ei in Ll, 

V • fir = 0 in D, 
fif = 0 on dQ, 

where e* denotes the ith standard canonical basis vector of M^. We also consider its finite 
element approximation {qi,, A/j) £ VhX that satisfies 


[ V(fif - fif;,) ; Vxdx - / (A - Xh)V ■ xdx = 0 G 14, 

J Q J Q 

/ caV • (fif - fir;,)dx = 0 Muj £ Mh. 


Using the definition of 5h we have 




\{dx^efi{z))i 


dxj dh^i ■ 6; 


dx 


and then 




Vg-.Ve^dx- [ W-eldx 
! JQ. 


^9h ■ - 


/ XhXI-eldx 

Jn 


Following the proof of Lemma we see that 


/ Vfif;, ; Vefidx - / A/,V • ef;dx = Ji + J 2 , 
In Jn 
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where 


Ji = ^ V{IhU-u) :Vgf^dx + j XhV ■ {hu - u)dx^ 

T&Th\T^ 

J 2 = -u + : Vgi/j dx+ A/iV • {hu -u + w'^)dx 


As in the proof of Theorem the estimates of Ji and J 2 follows from the properties of the 
correction functions and bounds for the Green’s functions that can be found for example in 
nziiin]. The estimate for the error of the pressure follows from similar arguments, we leave the 
details to the reader. □ 


Remark 2. Using the same arguments as Remark^ we can establish the same a priori error 
estimates given in Theorem^for u —{uh + wlf) and for p—(ph + w^). We note that the method 

interpolation operator, p^ + might not have average zero and a constant 0{h^~^^) must be 
added to make it in Lq(R). 


4.4) gives ph G M^, therefore, if Jh is the projection, then {p^ + wlf) G Lg(R), if is the 


5. Numerical Experiments 

We illustrate the performance of our method with some numerical examples. We consider the 
square domain R = [—1, 1]^, and we triangulate the domain with structured and non-structured 
triangular meshes. We tabulate the and semi-norm errors as well as the L°° and W^’°° 
semi-norm errors, with their respective order of convergence. Plots of approximate solutions 
and rate of convergence are also provided. For the case when the interface is not a straight line 
we need to integrate over curved region. We address this problem in Appendix [C] giving explicit 
quadrature formulas. 


5.1. Numerical examples for Poisson interface problem (1.1). Let u be the exact solution 
of problem (1.1), Uh be the solution by the method defined in (2.4). We define the error with 
respect to the Lagrange interpolant and the respective order of convergence (associated to 
the error and the norm) as follows 


C/i ■— Ufi IhU, 


\og{hi+i/hi) 


We will illustrate our results with two numerical examples using piecewise quadratic polyno¬ 
mials, i.e., k = 2. Note that Ih is then the piecewise quadratic Lagrange interpolant. 

Ex. 5.1.1 Consider an exact solution of problem 0 


u{x,y) 


f (2/3 + x) 3, ifx<l/3 
\ (4/3 —x)^, if X > 1/3. 


In this case, the interface T is a straight line. We summarize the results for structured 
and non-structured meshes in the following tables. 
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Ex. 


h 

e/i L2 

r 

Loo 

r 

l|Ve;,||L2 

r 

Ve/i LOO 

r 

3.5e-l 

8.41e-5 


1.83e-4 


1.53e-3 


3.72e-3 


1.8e-l 

7.49e-6 

3.49 

2.33e-5 

2.97 

2.76e-4 

2.47 

9.29e-4 

2.00 

8.8e-2 

6.63e-7 

3.50 

2.92e-6 

2.99 

4.92e-5 

2.49 

2.32e-4 

2.00 

4.4e-2 

5.85e-8 

3.50 

3.64e-7 

3.00 

8.74e-6 

2.49 

5.80e-5 

2.00 

2.2e-2 

5.16e-9 

3.50 

4.56e-8 

3.00 

1.55e-6 

2.50 

1.45e-5 

2.00 


Table 1. Errors and orders of convergence of method (2.4) for Ex. 5.1.1 


on 


structured meshes. 


h 

L2 

r 

e/i L°° 

r 

l|Ve;,||L2 

r 

Ve/i L°° 

r 

5.0e-l 

4.92e-4 


1.99e-3 


9.95e-3 


9.54e-2 


2.5e-l 

3.86e-5 

3.67 

2.37e-4 

3.07 

1.71e-3 

2.54 

2.41e-2 

1.99 

1.3e-l 

3.86e-6 

3.32 

4.42e-5 

2.43 

3.57e-4 

2.26 

7.98e-3 

1.59 

6.3e-2 

3.89e-7 

3.31 

5.99e-6 

2.88 

7.82e-5 

2.19 

1.92e-3 

2.06 

3.1e-2 

4.44e-8 

3.13 

7.74e-7 

2.95 

1.87e-5 

2.06 

5.06e-4 

1.92 


Table 2. Errors and orders of convergence of method (2.4) for Ex. 5.1.1 


on 


non-structured meshes. 


We observe optimal order of convergence for both, structured and non-structured 
meshes, confirming the results of Theorem and By structured mesh we mean 
bisecting in the northeast direction each cell of a uniform rectangular grid. When 
r is a straight line, the numerical quadrature is simplified, otherwise we need to 
integrate over curved elements (see Appendix]^. A superconvergence phenomenon 
is also observed in the norm for structured meshes, inherited from the problem 
without interface. 


).1.2 Consider an exact solution of problem 


u{x) = 


1 - log(3r), 


if r < 1/3 
if r > 1/3 


X G [-l,l]^ 


where r = ||x|| 2 . We summarize the errors and order of convergence in the following 
tables 


h 

L 2 

r 

e/i L°“ 

r 

l|Ve;,||L2 

r 

Ve/i L°° 

r 

7.1e-l 

1.70e-2 


4.76e-2 


1.89e-l 


4.37e-l 


3.5e-l 

1.73e-3 

3.29 

5.09e-3 

3.22 

3.66e-2 

2.37 

1.29e-l 

1.76 

1.8e-l 

1.49e-4 

3.54 

7.41e-4 

2.78 

5.66e-3 

2.69 

3.05e-2 

2.08 

8.8e-2 

1.22e-5 

3.61 

6.82e-5 

3.44 

8.48e-4 

2.74 

5.99e-3 

2.34 

4.4e-2 

1.16e-6 

3.40 

1.39e-5 

2.29 

1.54e-4 

2.46 

1.78e-3 

1.75 

2.2e-2 

1.09e-7 

3.41 

2.08e-6 

2.74 

2.71e-5 

2.51 

5.06e-4 

1.81 

l.le-2 

9.16e-9 

3.57 

2.39e-7 

3.12 

4.71e-6 

2.52 

1.36e-4 

1.90 


Table 3. Errors and orders of convergence of method (2.4) for Ex. 5.1.2 


on 


structured meshes. 


























HIGHER-ORDER FEM FOR ELLIPTIC PROBLEMS WITH INTERFACES 


19 


h 

L2 

r 

Loo 

r 

l|Ve;,||L2 

r 

Ve/i LOO 

r 

1.8e-l 

8.87e-5 


3.97e-4 


3.80e-3 


2.53e-2 


9.0e-2 

9.73e-6 

3.29 

7.46e-5 

2.49 

9.04e-4 

2.14 

7.43e-3 

1.82 

4.7e-2 

l.lle-6 

3.33 

1.06e-5 

3.00 

2.15e-4 

2.21 

2.58e-3 

1.63 

2.4e-2 

1.30e-7 

3.15 

1.42e-6 

2.95 

5.06e-5 

2.13 

7.34e-4 

1.84 

1.2e-2 

1.59e-8 

3.14 

2.24e-7 

2.76 

1.27e-5 

2.07 

2.16e-4 

1.83 

6.1e-3 

1.96e-9 

3.04 

3.15e-8 

2.85 

3.15e-6 

2.02 

5.55e-5 

1.98 


Table 4. Errors and orders of convergence of method (2.4) for Ex. 5.1.2 


on 


non-structured meshes. 



Figure 2. Plot of the L°° and W^’°° errors (left) and the approximate solution 
(right) by method (2.4) for Ex. 5.1.2, on non-structured meshes. 


In this example we are considering a curved interface, then quadrature rules over 
elements with one curved edge are needed. In Appendix we address this issue 
giving higher-order quadrature rules over this kind of elements. We observe optimal 
convergence, confirming our results in TheoremsandAs in the previous example, 
superconvergence phenomenon is observed in the case of structured meshes. 


5.2. Numerical examples for Stokes interface problem (1.2). In order to illustrate our 


method for the Stokes interface problem we consider the pair of inf-sup stable finite element 
spaces Vh and M^, piecewise quadratic polynomial for the velocity and piecewise constant 
polynomials for the pressure (P^ ~ Po)- We define the errors 


e“ := Uh - hu, el := ph - JhP 


With this finite element spaces, the value of /c is 1 (see assumption A2) , then is the 
piecewise linear Lagrange interpolant and Jh is the L?' projection on the spaces of piecewise 
constant polynomials. We note that the optimal order of convergence for velocity is 2, this is 
the reason we are using the piecewise linear Lagrange interpolant although we are using piecewise 
quadratic in our finite element space. 
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Ex. 5.2.1 Consider an exact solution of Stokes problem (1.2) on 11 = [—1,1] 


u{x,y) = 
p{x,y) = 


U2{x,y) = 


0 

U2{x,y) 

x'^ + y‘^ + 1/3, if X < 1/3, 
— 8/3, if X > 1/3 


2/3 + X, if X < 1/3, 
4/3 — X, if X > 1/3 

(x,y) e Q. 


In this case the interface is F = {(x,y) G 11 : x = 1/3}. 


h 

\K\\l^ 

r 

3.5e-l 

1.27e-2 


1.8e-l 

3.42e-3 

1.89 

8.8e-2 

8.84e-4 

1.95 

4.4e-2 

2.25e-4 

1.98 

2.2e-2 

5.66e-5 

1.99 


1/3'*^11 

r 


1.22e-2 


1.55e-l 

3.42e-3 

1.83 

7.86e-2 

9.67e-4 

1.82 

3.93e-2 

2.60e-4 

1.89 

1.96e-2 

6.84e-5 

1.93 

9.79e-3 


r 

l|Ve;}||L^ 

r 

0.98 

1.92e-l 

1.08e-l 

0.83 

1.00 

5.71e-2 

0.92 

1.00 

3.06e-2 

0.90 

1.00 

1.58e-2 

0.95 


Table 5. Errors and orders of convergence for velocity, on structured meshes. 


h 

II P II 

\\^h\\L‘^ 

r 

II P II 

r 

3.5e-l 

6.01e-2 


9.30e-2 


1.8e-l 

2.22e-2 

1.44 

6.07e-2 

0.62 

8.8e-2 

7.45e-3 

1.58 

3.52e-2 

0.79 

4.4e-2 

2.34e-3 

1.67 

1.92e-2 

0.88 

2.2e-2 

7.09e-4 

1.73 

l.OOe-2 

0.93 


Table 6. Errors and orders of convergence for the pressure, on structured meshes. 


We observe optimal convergence for the velocity and pressure, supporting our 
result in Theorem We also observe a superconvergence phenomenon for the 
error of the pressure. 


Ex. 5.2.2 Consider a exact solution of Stokes problem (1.2) on H = [—1,1] 


u{x,y) 


p{x,y) 


/ 

V 


ui{x,y) = I 
U2{x,y) = I 
4 - I if r < 


3y 

if r 

< 

1/3, 

3r y 

if r 

> 

1/3 

-3, 

if r 

< 

1/3, 

X 3r 

if r 

> 

1/3 

1/3, 





if r > 1/3 ’ 


\ 

/ 


for (x, y) G H and r = y^x^ + In this case the interface is the circumference of 
radius 3, i.e. F = {(x, y) G H : r = ^}. 
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h 

\K\\l^ 

r 

2.5e-01 

3.02e-2 


1.3e-01 

8.48e-3 

1.83 

6.3e-02 

2.03e-3 

2.06 

3.1e-02 

5.09e-4 

2.00 

1.6e-02 

1.26e-4 

2.02 


Wo'^W ^ 

r 


3.99e-2 


5.16e-l 

1.79e-2 

1.16 

2.79e-l 

5.35e-3 

1.74 

1.36e-l 

1.68e-3 

1.67 

6.84e-2 

4.22e-4 

1.99 

3.36e-2 


r 

II Ve^llioo 

r 

0.89 

8.10e-l 

5.48e-l 

0.56 

1.03 

3.35e-l 

0.71 

0.99 

2.06e-l 

0.70 

1.03 

1.05e-l 

0.97 


Table 7. Errors and orders of convergence for velocity, on structured meshes. 


h 

II P II 

l|e^llL2 

r 

II P II 

Wh\\L°° 

r 

2.5e-l 

1.39e-l 


1.84e-l 


1.3e-l 

3.39e-2 

2.04 

7.71e-2 

1.26 

6.3e-2 

1.32e-2 

1.36 

4.29e-2 

0.85 

3.1e-2 

3.79e-3 

1.80 

2.36e-2 

0.86 

1.6e-2 

1.46e-3 

1.38 

1.27e-2 

0.90 


Table 8. Errors and orders of convergence for pressure, on structured meshes. 



-1 -0.5 0 0.5 1 


Figure 3. Plot of velocity (left) and discontinuous pressure (right). 


We observe optimal convergence for the velocity and pressure, confirming our result 
in Theorem]^ We also observe from the plot of the approximate pressure (see Figure 
|Ex. 5.2.2| ) that there is not pollution beyond the interface. 

6. Final discussion and future work 


In this paper we have developed higher-order finite element methods for Poisson interface 
problem (1.1) and Stokes interface problem (1.2). We have proved optimal convergence results 
for the methods, recovering the convergence of the finite element method for a problem without 
interface. We have presented numerical experiments validating our theoretical results. The 
stiffness matrix remains the same as the problem without the interface. To the best of our 
knowledge, this is the first family of numerical methods that have provable higher-order accuracy 
(or arbitrary order) for these problems. 

The algorithms and theory developed point to a number of promising research directions 
and opportunities. In the future, we plan to use our techniques to develop and analyze high 
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contrast problems for elliptic problems with discontinuous diffusion coefficients. Additionally, it 
is natural to move in the direction of evolution problems, such as the time-dependent version of 


problem (1.2) with moving interfaces. We plan to develop time stepping methods where in each 


time step we used the methodology developed here to solve the resulting problem. Navier-Stokes 
equations and more challenging problems such as fluid-structure interactions will be explored in 
the future. 

Finally, more straightforward generalizations are possible. In this paper, for convenience, 
we used triangular meshes, but it is clear that rectangular or quadrilateral meshes can be 
used. Indeed, the correction function did not use the fact that the mesh is triangular. Three 
dimensional problems should not pose a serious difficulty. The key to our method is the correction 
functions w^. The construction of corresponding correction functions in three dimensions seems 
natural. Details will appear in the Ph.D. thesis of the second author. 


Appendix A. Proof of estimate (3.17) 


We first use a dyadic decomposition of D. D = D* U Uj=i where Qj := {x G D : dj+i < 
|x —zl < dj}, and dj = 2~K Here, J = log 2 (I//i) and H* = Vt.r\Bh where Bh is the ball of radius 
h centered z. 

Let us use the notation Dj = Dj n dRh\dVt and D* = n dRh\dQ then we have 

,7 

\\Dng\\L^{dRh\dQ) = \\Dng\\L'^(D*) + 

J = 1 

First we bound the first term. Using Cauchy-Schwarz inequality and the fact the one dimen¬ 
sional measure of D* is h we have 

\\Dng\\L^{D*) ^ Ch^^^\\Vg\\L2i^D*y 
If we use a trace inequality then we have 


h^^‘^\\'^g\\D{D*) < + h\\D‘^g\\D{n*nBh)) < C{\\Vg\\L2(^Q-j -|- h\\D‘^g\\L2(^Q-j). 

Elliptic regularity will give 

< C!{\\6h\\L^(n) + h\\6h\\H^(n)) < 

which shows that 

C 

11 -^^^ 5111 . 1 ( 0 *) < 

To take care of the sum we use Cauchy-Schwarz inequality again and get 


J J 


7=1 


7=1 


Using the Green’s function representation of g and the fact that Dj is distance dj away from 
z one can show (see m) 

C 

W'^gWi^iDj) < -J-- 


Yj \\^r^g\\L^Dj) < 

7=1 7=1 ^ 


7=1 


0(2-^+^ 


1 )< 


c 

~h' 


Hence 
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Appendix B. Derivation of jumps (4.2) 

We define for e > 0 the set = {x G D : dist{x,T) < e}, where dist is the standard 
distance function. We will also denote From equation (4.1a), taking divergence, 

multiplying by an arbitrary G of compact support in and integrating we obtain 

/ V • Vp(j)dx = V • f4>dx + / V • B(j)dx 

JUf: JUf: 

(B.l) = [ V-f^dx- [ f3-V(j)ds. 

Jn^ Jr 

For the left-hand side of the equation above we integrate by parts twice obtaining 


/ V • Vp 4>dx = / Dp'^ ■ n(j)ds -|- / Dp ■ n(j)ds 

In. Jr+ Jrr 


— / D(f)-np'^ds— / Dcj) ■ np ds+ pV ■'S/4>dx. 

Jr+ Jrz Jn^ 

By the definitions of F^ we can fix its normal vectors by the normal vector to the interface F. 
Then, taking the limit as e goes to zero we obtain 


(B.2) 

lim 

/ pV 

o 

II 

> 

e^O ^ 

In^ 


(B.3) 

lim 

[ 

< 

II 

e^O ^ 

In^ 



I [Dp ■ n] 4>ds — J Dcj) ■ n\p]ds 

We need to write the right-hand side in ( B.l| ) in terms of the normal and tangential derivative. 
Let 9 be the angle between the x direction (x-axis) and n direction and let R{9) the rotation 
matrix defined by 

cos(6*) — sin(0) 

sin(0) cos(0) 


R{9) = 


Then, V(/> = (Dncj), Dt(l>)R{9Y. Using this equivalency in (B.l) and integration by parts we 
obtain 

[ P • Vcj)ds = [ p- {(Dr^cj), Dt(t>)R{ef) ds 


= J^{pR{e))-{Dr,cj),DtJ))ds 


= J P ' iDn4>, Dt(i))ds 

= [ PiDn(j)ds - [ ^Mds 

Jr Jr ds 


Here we have defined p = PR{9). Matching this result with the limits in (B.2) we can conclude 
that 


(B.4) 


b] = h 


[o„Pl = 


Now, to get the normal jump of the gradient of the velocity we applied equation (1.2d), 
resulting 


(B.5) 


[Dnu] = P + [pn] = P + Pin. 
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Appendix C. Quadrature over curved elements. 

The method defined in (2.4) requires the integration, over an element T G of • Vf. 
Now, since the correction function is defined as piecewise polynomial of degree k on each 
piece of the restriction = T n we need some quadrature formulas on this region, say 
curved elements (polygons with one edge curved). To do so, we observe that Lt divides T in 


two regular polygons T^, one is a triangle and the other a quadrilateral, see figure]^ In order 
to make the presentation clearer, we will write the integral in terms of the coordinates system 
(r, rj) (associated to the vector r and r/, with ij pointing outward T~) such that the origin is 
the midpoint of Lt, the point . Let / G then we write 


(C.l) 


/t± 


f{x,y)dxdy = 


r /■'■t/2 niT) 

/ f{x,y)dxdy± / / f{T,r])dr]dT 

lf± J-ttIi jo 


where, 7 (r) is the equation for the curve T on T. The first integral in (C.l) is over a polygon 


(triangle or quadrilateral) and we can use quadrature rules over polygons to approximate it, 
whereas the second integral is related to the curve, and we will use some Gaussian quadratures. 
First step, we take n Gaussian points over the segment Lt, equivalently, the interval 

[—r'r/ 2 , rr/ 2 ] in the (r, 7 /)-axis, and for each point we compute its projection onto the curve T 
in the 77 —direction, say 7 (rj). Now we compute, the Gaussian points for the segment [ 0 , 7 (rj)], 
we denote them by {i]i,j}'j=i- Therefore 


rrT/2 nir) 

Jo 


fl \/l ^ I7(t)| ^ 

f[T,y)dTdr] ^ 

i=l 




where iOi are the weights in the interval [—1,1]. See hgurej^for illustration of this notation. See 
also figure 



Figure 4. Illustration of the Gauss points. 
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